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(S| Abstract 



A technique using vector Slepian harmonics and vector Slepian multipole fields is presented for a general treatment of the inverse 
(S| problem of high numerical aperture focusing. A prescribed intensity distribution or electric field distribution in the focal volume 
Q^is approximated using numerical optimization and the corresponding illuminating field at the entrance pupil is constructed. Three 
examples from the recent literature are chosen to illustrate the method. 
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1. Introduction 



. ^ Engineering the focal spot of a high numerical aperture (NA) 
^ ^'lens is a central problem in several optical applications such as 
Q laser scanning microscopy, optical trapping and laser microma- 
^ chining. To obtain the desired spot, one usually has to solve a 
Q jiumerical optimization problem. There exist several examples 
* ^ for the design of annular amplitude or phase masks that tailor a 
given input beam for a specific application [1-4]. These calcu- 
lations usually rely on the Debye-Wolf diff'raction integral [5] 
, ^^ to determine the focused field. 

Instead of the direct numerical evaluation of the Debye-Wolf 
CO integral, some methods use analytical series expansions. These 
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methods include the multipole theory of focusing [6], the tech- 
nique of Kant using Gegenbauer polynomials and spherical 



"Bessel functions [7], the extended Nijboer-Zemike (ENZ) ap- 
fNj proach [8] or the scalar eigenf unction expansion of Sherif et al. 
using Slepian's prolate spheroidal functions [9]. Some of these 
methods were also used for solving the inverse problem of fo- 
CS| cusing [10-12]. Recently, a technique using dipole arrays has 

also been proposed [13, 14]. 
^ The methods above, however, possess some inherent draw- 
backs. The scalar Slepian method [9] is not suitable for opti- 
mization problems with a prescribed intensity pattern [12], the 
dipole array method is tailored for specific problems [13, 14] 
" and the multipole theory [6] lacks directionality thus requires 
additional constraints to ensure directional illumination. Al- 
though the ENZ theory [8] is suitable for aberration retrieval 
where a linearized intensity approximation can be used [11], it 
becomes, similarly to the method of Kant [7], computationally 
challenging for arbitrary focal intensity profiles. 

Based on the multipole theory of focusing [6], we have re- 
cently proposed an orthonormal basis of vector Slepian har- 
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monies (VSLHs) which are naturally suitable for approximat- 
ing the illumination in a high NA system [15]. Its main ad- 
vantage is that a subset exhibits excellent directionality, i.e. its 
angular energy distribution is confined to the solid angle of il- 
lumination. Each VSLH basis function represents the vectorial 
plane- wave amplitudes (PWAs) of a focused field described by 
a corresponding vector Slepian multipole field (VSLMF). The 
directionality of the VSLHs allows us to approximate common 
illumination PWAs using a smaller number of coefficients than 
a representation using the vector spherical harmonics does [15]. 

In this paper, we demonstrate the applicability of these novel 
vector bases for general inverse problems in high NA focusing. 
As illustrations, three examples are considered: the optical nee- 
dle [13], tube [14] and bubble [4]. 



2. The theory of vector Slepian harmonics and multipole 
fields 

First we give a brief introduction to the new VSLHs and 
VSLMFs discussed more comprehensively in Ref. [15]. One 
can approximate a divergence-free, monochromatic focused 
electric field that satisfies the vector Helmholtz equation as 



Tnax 

E(r) = ^^c^,-0^,-(r), 



(1) 



where the vector functions 0^/(r) are the VSLMFs and Cmi de- 
notes complex expansion coefficients (determined in this pa- 
per by numerical optimization, as discussed in Section 3). The 
VSLMFs themselves are linear combinations of a finite number 
of vector multipole fields (VMFs) M/;^(r) and N/;^(r) [15] of the 
same order m, i.e. 
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„<r) := ^ [uiiJ^imir) + vaMm(r)] , (2) 



l^im 



where L denotes the maximal degree of the contributing VMFs 
and the summation starts at 

^ ifm = 0, 

^ \\m\ otherwise. 

On the calculation of the coefficients uum and vum of Eq. (2), see 
Ref. [15]. 

The azimuthal properties of 0^/(r) are solely determined by 
the order m. In most cases, e.g. for focused fields of linearly 
or circularly polarized input beams or cylindrical vector beams, 
only VSLMFs of certain values of m contribute to the focused 
field. The second integer / indexes VSLMFs of the same order 
m, however, the limit /max of Eq. (1) is usually chosen to be 
much smaller than the number A^^ of available VSLMFs for a 
given order m, where 



Nm := 



\L ifm = 0, 

L - Iml + 1 otherwise. 



(3) 



To clarify the reason for this choice, we must first recall the 
plane- wave representation of these fields. 

The VSLMFs can also be expressed as a superposition of 
elementary plane- waves, i.e. 
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where k = In I A denotes the wavenumber {A is the wave- 
length), 6s and (ps are angular spherical coordinates, and § = 
sin 6s cos (ps + sin 6s sin (ps % -h cos 6s is the propagation unit 
vector of a plane- wave component with an amplitude and polar- 
ization specified by Vmiifis, (ps)- Between JJ mi(6s, (ps) and s the 
trans versality relation \]mi(6s, (ps) = holds. 

Analogously to Eq. (2), the functions \Jmi(6, (p) are also de- 
fined as a linear combination: 

L 

lJnii(6, (p) := ^ [UilmYlm(6, (p) + VnmZimifi, 0)] , (5) 

where Yim(6,(p) and Zim(6,(p) are vector spherical harmonics 
(VSHs) [15]. Since (i) the VSHs represent the PWAs of the 
VMFs [16] similarly to Eq. (4) and (ii) VSLMFs are defined 
as a linear combination of the VMFs according to Eq. (2), the 
same expansion coefficients uum and vum appear in Eq. (5) as 
in Eq. (2). Furthermore, following from the linearity of our 
expressions, the combination 



'■max 



A(6s, (Ps) = 2_jZj Cmi^miiOs, (ps) 
m i-l 



(6) 



of VSHs with the coefficients Cmi of Eq. (1) gives the PWAs 
A(^^, (ps) of the focused field. The electric field at the entrance 
pupil of the lens can then be found using the vectorial ray- 
tracing method of Torok et al. [17]. 

Although Eq. (4) resembles the Debye-Wolf integral [5] 
closely, it is important to stress that the integration is performed 



over the unit sphere of all possible plane-wave directions, not 
just for the spherical cap 

:= (ps)\0<6s< sin-^NA), < 0, < In}. (7) 

To be applicable to systems whose NA restricts possible plane- 
wave directions to Sq, the VSLHs l]mi(6, (p) are constructed to 
maximize T]mi, the fraction of the energy falling onto Sq [15]. 
The VSLHs \Jmi(6, (p) for a given m are ordered by decreasing 
rimi, thus those with lower / have less "energy leakage" outside 
Sq. In our examples, we restrict ourselves to basis functions 
with ^m/ ^ 98%. This leads to choosing the maximal index /max 
to be less than A^^^. We stress the fact that using this method, 
no extra constraints are necessary to ensure the directionality 
of the illumination, which would be inevitable when using the 
multipole theory of focusing [6]. 



3. The inverse problem 

Next we demonstrate how to use our method of VSLMFs and 
VSLHs to treat the inverse problem, i.e. to find the illumination 
for a prescribed focal intensity profile. Here we only consider 
examples of cylindrical vector beams that can be described as 
a superposition of radially and azimuthally polarized fields. In 
the context of VSLHs this means a choice of m = [15]. How- 
ever, all optimization tasks could be performed in a straight- 
forward way for linear, circular or more general polarization 
states as well, merely with different choices of m. Naturally, 
our method is easily applicable to inverse problems on the elec- 
tric field distribution, too, provided that a physically feasible 
focal electric field is prescribed. 

In our calculations an aplanatic lens with NA = 0.95 was 
used. Assuming a cylindrical coordinate system (p, (p, z) cen- 
tered at the focus, a rotationally symmetric intensity distribu- 
tion /(p, z) was prescribed in the focal volume. 

The focused field can be written as 



'■max 

E(r;Q) = ^QOo.-(r), 



(8) 



where the task is to find the coefficients c/ that minimize the 
cost function 



F{ci) := 



[r(r;Q)-E(r;Q)-/(r)]' d^r. 



(9) 



the squared error of the intensity integrated over some reason- 
ably chosen volume V around the focus (the asterisk denotes 
the complex conjugate and we have simplified the notation by 
omitting the zero m index from cq/). 

The optimization was performed by sequential least squares 
programming [18]. Random initial values of c/ were chosen and 
the process converged to a local minimum of F(ci). For a given 
problem, this optimization process was repeated 40 times and 
the lowest minimum was selected. Although not guaranteed to 
find a global optimum, the process yielded satisfactory results 
from a practical point of view. If attaining the global optimum 
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is of crucial importance, computationally intensive global opti- 
mizing algorithms can be used. We note that 40 runs generally 
yielded only at most 3 different local minima and there was less 
than 2% difference in the value of the cost function for the dif- 
ferent local minima. 

For the examples presented in this paper, we only considered 
real values of c/ that correspond to a binary phase modulation of 
the radial and azimuthal components of the input beam. When 
allowing complex- valued c/ (i.e. a locally elliptically polarized 
input beam), the size of real optimization parameters doubled 
and the computation time increased by a factor of 2.4-2.7 on 
the average. However, we did not obtain any better local optima 
in terms of cost function value for the examples considered. 

The PWAs of E(r) can be obtained as 



■•max 



(10) 



i=l 



Since rjmi < 1 [15], a small amount of energy leakage may still 
occur outside Sq- Because of that, the exact focused electric 
fields were also calculated by Debye-Wolf integration using 
Chirp Z-Transform [19] after enforcing a hard limit at the edge 
of the entrance pupil. 

The prescribed focal intensity functions for the needle, tube 
and bubble were 



(11a) 
(lib) 
(11c) 



I(p) oc exp|-3-), 
/(p)oc- exp^-- 
I(p, z) oc — exp 



p' + (z/3)2 



f,2 



respectively, where the parameters a and b were chosen such 
that the full width at half maximum (FWHM) of all func- 
tions in the focal plane was OA A. For the optical tube and 
bubble this value refers to the inner FWHM size of the dark 
spot. For the optical bubble, a physically feasible aspect ra- 
tio FWHMJFWHMp = 3 was chosen, as seen in Eq. (11c). 
The integration in Eq. (9) was performed over a cylinder with 
a length of 10 A and a radius of 2 A for the needle and the tube, 
and over a prolate ellipsoid of revolution with a length of 12 A 
and a radius of 2 A for the bubble. The shape of the integration 
volumes was chosen to be similar to each prescribed focal spot; 
and their lateral size was large enough so that the prescribed in- 
tensity at the edges was at least 9 orders of magnitude smaller 
than the peak value. 

In the case of the optical needle and tube, the axial length 
of the integration volume approximately matched the maximal 
length of the focal spot attainable by the VSLMF basis func- 
tions involved in the process. This maximal length can be in- 
creased by increasing L of Eq. (2) [15]. For all examples pre- 
sented here, we chose L = 30. 

For both the optical needle and tube, the 1 1 most concen- 
trated Oo/(r) functions were used. For the needle we included 
only Oo/(r) corresponding to radially polarized illumination. 



because the azimuthally polarized components would only in- 
crease the lateral FWHM and would not contribute to the on- 
axis electric field. For analogous reasons, only Oo/(r) corre- 
sponding to azimuthally polarized illumination were included 
for the optical tube. The bubble included both types of polar- 
izations. Moreover, since the ideal optical bubble demands zero 
electric field in its center, we introduced the extra constraint 
E(r = 0; c/) = in the optimization process. 

4. Results 
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Figure 1: Radial and axial intensity cross-sections of the optical (a)-(b) nee- 
dle, (c)-(d) tube and (e)-(f) bubble, respectively, with the solid line showing 
the result of the Debye-Wolf integration and the dashed line the VSLMF ap- 
proximation (everywhere except in (b) the two curves practically overlap). The 
insets of show intensity contour maps in (a), (c), (e) the focal plane and (b), (d), 
(f) the meridional plane, respectively. 

Fig. 1 depicts the intensity distributions obtained from the 
optimization process for each spot and Fig. 2 shows the corre- 
sponding electric fields at the entrance pupil of the lens. In the 
case of the needle (Figs. l(a),(b)), the length of the spot is 9.09/1 
(axial full width at 90% maximum), while the lateral FWHM is 
0.41/1 in the focal plane and below 0.44/1 throughout its full 
length. The beam purity (expressing the fraction of energy of 
the longitudinally polarized component inside the needle vol- 
ume, defined in Ref. [13]) stays above 78% along the needle 
with a value of 86% in the focal plane. Hence, we can con- 
clude that the quality of the needle resulting from our method is 
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Figure 2: The intensity distribution and polarization of the electric field at the entrance pupil for the optical (a) needle, (b) tube, and (c) bubble {R is the radius of 
the entrance pupil). 



comparable to that of Wang et al. [13]. Fig. 1(b) shows a slight 
difference between the Slepian approximation of the focused 
field using the VSLMFs Oo/(r) and the result of the Debye- 
Wolf integration. This discrepancy can be attributed to the fact 
that 0.54% of the total energy falls outside the spherical cap of 
illumination. 

The optical tube shown in Figs. 1(c), (d) has a length of 9.09/1, 
measured along p = Pmax, where pmax is the peak radius in 
the focal plane. The lateral inner FWHM is 0.46/1 in the fo- 
cal plane and stays below 0.50/1 throughout the entire length of 
the tube. Again, our results are comparable with those of Wang 
et al. [14]. Here the energy leakage, when using the Slepian 
approximation, is only 0.04%, since the illumination intensity 
is low near the edge of the entrance pupil (as seen in Fig. 2(b)). 

We note that by increasing L the number of total VSLHs in- 
creases, including the number of those with rimi ~ 1- While 
more optimization parameters are introduced, the maximal at- 
tainable length of the tube and needle also increases, as already 
noted in the previous section. Thus when increased depth of 
field is of primary importance, a larger value for L can simply 
be chosen. 

Finally, the optical bubble of Figs. 1(e), (f) qualitatively re- 
produces the result of Bokor and Davidson [4]. In this case, the 
illumination (Fig. 2(c)) is constructed from 22 basis functions. 
The axial and lateral FWHM of the dark core of the bubble 
are 1.22/1 and 0.43/1, respectively. As seen in Figs. 1(e), (f), the 
intensity is zero at the center as prescribed, and the axial and 
transverse peak intensities are nearly equal in magnitude, which 
is advantageous e.g. in fluorescence depletion microscopy [4]. 
Again, the illumination intensity is low near the edge of the en- 
trance pupil, leading to an energy leakage of only 0.03%. 

5. Concluding remarks 

We have demonstrated with three well-known examples that 
VSLHs and VSLMFs are highly suitable for inverse problems 
of high NA focusing. Knowing the desired 3D intensity distri- 



bution in the focal volume, our method can be used efl'ectively 
to design the electric field distribution at the entrance pupil. 

Although the realization of the input field requires careful 
control of amplitude, phase and polarization, some practical 
methods using liquid crystal spatial light modulators or form- 
birefringent spatially variant subwavelength gratings have al- 
ready been proposed to achieve this task [20-22]. 

Finally, it is important to note again that VSLHs and 
VSLMFs are suitable for treating the inverse problem not only 
for the focal intensity but for the focal electric field as well. In 
fact, in that case a least squares approach similar to Eq. (9) leads 
to an easily solvable system of linear equations. 
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